#!/usr/bin/python

# a lil' sw to filter spectra if they have a 2 Da doublet...

import xlh
from random import Random
import sys, getopt


def exitOnError(errmsg, exit):
	sys.stderr.write('The following error(s) occurred:\n')
	sys.stderr.write(errmsg + '\n')
	sys.exit(exit)
def usage():
	print sys.argv[0], """-- a tool to validate OMSSA results in a X-link framework
	
Usage:

xfilter.py [-h] [-r RATIO] [-t TOL] [-u Unit] [-o File] [-d doublet] [-n INTF] -i File.mgf

 -h:	Print this help and exit
 -i:	Input file
 -r:	Peak ratio [Default: 0.85]
 -t:	Tolerance for ion matching [Default 0.5 Da]
 -u:	Tolerance unit [Default: Da]
 -o:	Output file name [Default: StdOut]
 -n:	Intensity Filter [Default: 3]
 -d:	Doublet distance [Default: 2]
	"""
	sys.exit(1)

Doublet = 2.0
Tolerance = 0.5
TolUnit = 'Da'
intRatio = 0.85
intFilter = 0.03
Fileout = '_fsp' + str(Random().getrandbits(64)) + '.mgf'
Filein = ''

try:
	optlist, args = getopt.getopt(sys.argv[1:], 'hi:r:t:u:o:n:d:')
except getopt.GetoptError:
	usage()

if not optlist:
	usage()

for o, a in optlist:
	if o == '-h':
		usage()
	if o == '-r':
		intRatio = float(a)
	if o == '-d':
		Doublet = float(a)
	if o == '-n':
		intFilter = float(a) / 100
	if o == '-t':
		Tolerance = float(a)
	if o == '-u':
		TolUnit = a
		if TolUnit not in ['Da', 'ppm']:
			exitOnError(TolUnit + " is not a valid unit, asshole!", 1)
	if o == '-o':
		Fileout = a
	if o == '-i':
		Filein = a
		
		
# do some checks
if not Filein:
	exitOnError('Missing input file', 1)
	
	
# some other in the future....
#### End of option parsing section

# open a mgf file:

spFile = xlh.MSAnalysis()
towrite = xlh.MSAnalysis()

spFile.parseMgfFile(Filein)

for charge in [3,4]:
	for sp in spFile.chargeSearch(charge).items():
		spsl = sp.spectrumSlice(min = sp.precursorMz()).filterOnIntensity(sp.maxIntensity() * intFilter)
		for peak in spsl.items():

			if TolUnit == 'Da':
				rtol = Tolerance
			elif TolUnit == 'ppm':
				rtol = peak[0] / 1000000 * Tolerance

			clpeaks = filter(lambda x: x[0] - peak[0] <= (Doublet + rtol) and x[0] - peak[0] >= (Doublet - rtol), sp.peaks()) 

			if filter(lambda x: x[1] / peak[1] >= intRatio or peak[1] / x[1] <= 1 / intRatio, clpeaks):
				towrite._spectraList.append(sp)
				break
			
			# if a doublet in 2 Da has been found stop here, otherwise check if there's a doublet in 4 Da...
			
			clpeaks = filter(lambda x: x[0] - peak[0] <= (2*Doublet + rtol) and x[0] - peak[0] >= (2*Doublet - rtol), sp.peaks()) 

			if filter(lambda x: x[1] / peak[1] >= intRatio or peak[1] / x[1] <= 1 / intRatio, clpeaks):
				towrite._spectraList.append(sp)
				break
			
towrite.dumpMgfFile(Fileout)



